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Abstract 

The use of Mean-Field theory to unwrap principal phase patterns has been recently 
proposed. In this paper we generalize the Mean-Field approach to process phase 
patterns with arbitrary degree of undersampling. The phase unwrapping problem 
is formulated as that of finding the ground state of a locally constrained, finite size, 
spin-L Ising model under a non-uniform magnetic field. The optimization problem 
is solved by the Mean-Field Annealing technique. Synthetic experiments show the 
effectiveness of the proposed algorithm. 
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1 Introduction 



The determination of the absolute phase from a fringe pattern is an important 
problem that finds applications in many areas: homomorphic signal process- 
ing [1], solid state physics [2], holographic interferometry [3], adaptive or com- 
pensated optics [4], magnetic resonance imaging [5] and synthetic aperture 
radar interferometry [6]. 

In all these applications one obtains a two-dimensional fringe pattern whose 
spatially-varying phase is related to the physical quantity to be measured. The 
computation of phase by any inverse trigonometric function (e.g. arctangent) 
provides only principal phase values, which lie between ±7r radians. The pro- 
cess of phase unwrapping (PU), i.e. the addition of a proper integer multiple 
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of 2tt to all the pixels, must be carried out before the physical quantity can 
be reconstructed from the phase distributions. Since many possible absolute 
phase fields are compatible with a given fringe pattern, phase unwrapping is 
ill-posed in a mathematical sense: Hadamard [7] defined a mathematical prob- 
lem to be well-posed if a unique solution exists that depends continuously on 
the data; in this case, the uniqueness requirement is violated. 

Ill-posed problems arise frequently in many areas of science and engineer- 
ing; well-known examples are analytic continuation, the Cauchy problem for 
differential equations, computer tomography, and many problems in image 
processing and machine vision that involve the reconstruction of images from 
noisy data. 

The fact that a problem is not well-posed does not mean that it cannot be 
solved: rather, in order to be solved it must be first regularized by introducing 
additional constraints (prior knowledge) about the behaviour of the solution. 
Variational regularization corresponds to modeling the physical constraints 
of the problem by a suitable functional; the solution is then sought as the 
minimizer of this functional. 

In the last years, an increasing interest has been devoted to adapt methods 
from Statistical Mechanics to nonconvex optimization problems arising from 
the variational regularization of ill-posed problems. Geman and Geman [8] 
suggested that the Ising model is applicable to image restoration through the 
Bayesian formalism. This problem corresponds to searching the ground state 
of a finite-size Ising model under a non-uniform external field. Geman and 
Geman applied this to the recovery of corrupted images by using simulated 
annealing of a spin-S Ising model. After that, Gidas [9] proposed a new method 
based on a combination of the renormalization group technique and the simu- 
lated annealing procedure; then, Zhang [10] introduced Mean-Field Annealing 
to treat the image reconstruction problem, while Tanaka and Morita [11,12] 
applied the cluster variation method. Methods of statistical mechanics have 
also been used to study combinatorial optimization problems (see, e.g., [13] 
and references therein). 

In a couple of recent papers, the phase unwrapping problem was handled by 
methods of Statistical Mechanics. Simulated annealing was applied in [14]. 
In [15] the problem is solved by the Mean-Field Annealing (MFA) technique; 
PU is formulated as a constrained optimization problem for the field of integer 
corrections to be added to the wrapped phase gradient in order to recover the 
true phase gradient, with the cost function consisting of second order differ- 
ences, and measuring the smoothness of the reconstructed phase field. This is 
equivalent [15] to finding the ground state of a locally-constrained ferromag- 
netic spin-1 Ising model under a non- uniform external field. The optimization 
problem is then solved by MFA and consistent solutions are found in difficult 
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situations resulting from noise and undersampling. 

Mean Field Annealing is closely related to Simulated Annealing [8]. Both ap- 
proaches formulate the optimization problem in terms of minimizing a cost 
function and defining a corresponding Gibbs distribution. Simulated Anneal- 
ing then proceeds by sampling the Gibbs probability distribution as the tem- 
perature is reduced to zero, whereas MFA attempts to track an approximation 
of the mean of the same distribution. 

The algorithm described in [15] was constructed under the assumption that 
the possible values for the correction field were restricted to belong to the set 
{ — 1,0,1}. In this paper we generalize the MFA algorithm to the case of a 
set {— L, ... , L}, with L an arbitrary integer. The corresponding statistical 
system is a locally-constrained spin-L Ising model. 

The present generalization allows the processing of input phase patterns with 
arbitrary degree of undersampling; our experiments on synthetic phase fields 
show the effectiveness of the proposed algorithm. 

The paper is organized as follows. In sect. 2 the phase unwrapping problem 
is introduced and its ill-posedness is highlighted. In sect. 3 the deterministic 
MFA algorithm is described in detail. Then, in sect. 4 some experimental 
results on simulated phase fields are presented. Some conclusions are then 
drawn in sect. 5. 



2 Phase Unwrapping Problem 

We briefly recall here the phase unwrapping terminology, and refer the reader 
to [16] for a complete discussion. 

Given an absolute phase pattern f(x,y) on a two-dimensional square grid, 
what is actually measured is the wrapped phase field g(x,y) which can be 
expressed in terms of the / field through a wrapping operator, Wr, defined so 
that g(x,y) always lies in the interval [— it, +7r): 

g(x, y) = Wr[f(x, y)] = arg {exp[i/(x, y)}} . (1) 

Phase unwrapping means recovering the absolute phase field /, which is usu- 
ally related to the physical quantity to be measured, from the knowledge of 
the g field. This can be done in practice by estimating the absolute phase 
gradient from the wrapped phase field and integrating it throughout the 2-D 
grid. This simple method is effective only in absence of phase aliasing, i.e. if 
the phase field is correctly sampled. 
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In fact, if the Nyquist condition: 



|V/(x,y)|<7r, (2) 

where V is the discrete gradient, is verified everywhere on the grid, the ab- 
solute phase gradient is obtained by wrapping the gradient of the wrapped 
phase field, according to the formula: 

A(x,y)=Wr[Vg(x,y)]. (3) 

As mentioned, if condition (2) is satisfied, one has: 

Vf(x,y)=A(x,y), (4) 

and the /-pattern is obtained by integrating A along any path connecting all 
sites on the grid. 

The Nyquist condition is often violated because of undersampling of the sig- 
nal from which the principal phase is extracted. This can result either from 
system noise, or from critical values of the slopes of the physical surface which 
is analyzed through interferometry. For example, in the case of SAR interfer- 
ometry, the surface is the portion of Earth imaged from the sensor (usually 
air- or satellite-borne), while noise can arise from sensor thermal electronic 
motion, or from other sources of electronic signal disturbances. 

If the Nyquist condition is not satisfied everywhere on the grid, then the 
wrapped gradient A of the wrapped phase field is not assured to equal the 
absolute phase gradient. In this case, a more general relation must be written, 
rather than (4), i.e.: 

Vf(x,y) = A(x,y) + 27rk(x,y), (5) 

where k(x, y) is a vector field of integers. In this case, solving PU amounts to 
finding the correct field k. 

Phase aliasing conditions imply that the integration of field A depends on 
the path. The sources of this nonconservative behaviour are detectable by 
calculating the integral of the field A over every minimum closed path, i.e. 
the 2x2 square having the site (x,y) as a corner: 

y) = [A x (x, y) + A(y(x + 1, y) - A x (x, y + 1) - A y (x, y)] . (6) 

One can show that the integral I(x,y) will always have a value in the set 
{ — 1, 0, 1}. Locations with 1^0 are called "residues". In presence of residues, 
the field A is no more irrotational; this causes the path-dependence of the 
integration step previously described. 
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To restore the consistency of the phase gradient, then, the k field must satisfy 
the following consistency condition: 



Vx [A(x,y) + 2nk(x,y)] = 0, (7) 

where V x • is the discrete curl operator. Since there are many possible k 
fields satisfying eq. (7), PU is an ill-posed problem according to Hadamard's 
definition. 

One of the most classical and widely-used algorithms for phase unwrapping 
is the Least Mean Squares (LMS) approach [17], which consists in finding the 
scalar field / whose gradient is closer to A in the Least Squares sense, i.e. the 
minimizer of: 

E(W-A) 2 . (8) 

As mentioned before, in [15] a variational approach has been used, and the 
field k was "chosen" as the minimizer of the following functional: 

R = ^£[Vx/(z + l,y)-V*/(:r,y)] 2 

+ 4^ E [ v */(^ V + 1 )- v -/(^ y)f 

+^ E [V s /(x +1,2/)- V y f(x, y)f , (9) 

subject to constraint (7). Due to (5), R is a functional of k, i.e. R = R\k], and it 
measures the smoothness of the reconstructed phase surface. The optimization 
problem was then solved by a MFA algorithm under the assumption that the k 
field be restricted to take values in { — 1, 0, 1}. In the next section we generalize 
the MFA algorithm to the case of k fields belonging to {— L, ... , L}, with L 
an arbitrary integer. 



3 The algorithm 



As explained in sect. 2, we assume the solution of PU to be the minimizer of the 
functional (9), subject to constraint (7). Let us assume that the possible values 
of the k field are restricted to belong to {— L, ... , L}. The field k may then 
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be regarded as a system of spin-L units. We assume the Gibbs distribution: 

P[k] 



exp 




T 




Ek' exp 


T 



(10) 



where the sum is over the k' fields satisfying (7), and T is the statistical 
temperature. Inconsistent fields k' are assumed to have zero probability. 

Following Mean-Field theory [18], we consider a probability distribution for 
the correction field k which treats all the variables as independent, i.e. it is 
the product of the marginal distributions of each variable. 

Let p x (x,y,a) be the probability that k x (x,y) = a, with a = —L,... ,L, 
and p y (x, y, a) be the corresponding probability for k y (x, y). Normalization of 
these marginal probabilities implies a penalty functional: 



©W = E 

(x,y) 



V x (x, y) ^1 - Px(x, y, a)j + V y (x, y) ^1 - ^ p y (x, y, a) j 



11^ 



where {V} are Lagrange multipliers. 

The entropy of the system, in the mean field approximation, is: 

L 

%] = - E H [Px(x, y, a) log p x (x, y, a) + p y (x,y, a) log p y (x,y, a)]. (12) 

(x,y) a=-L 



It is useful to introduce the average fields m = (k) p and Q = (k 2 ) p , defined 
by: 



m x (x,y)= apx(x,y,a), m y (x,y) = a Pyi x ^y^ a )\ ( 13 ) 

a=— L a=—L 

L L 

Qx(x,y) = E a 2 px(x,y,a), Q y (x,y) = ^ a 2 p y (x,y,a). (14) 



a=— L 



a=—L 



The average U of the cost functional R is called internal energy. It is easy to 
show that the internal energy depends only on m and Q: 



U[m, Q] = (i?[A + 27rk]) p . 
A penalty functional is introduced to enforce constraints (7): 



(15) 
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r ( m ] = Yl v) [ m x( x i v) + m y( x + !> y) 

-m x (x,y + 1) - m y (x,y) + I(x,y)\, (16) 
where {A} is another set of Lagrange multipliers. 

Let us now introduce an effective cost functional, the free energy, which de- 
pends on T: 

F[p] = U[m, Q] - TS[p] + r[m] + 9[p] (17) 



The free energy is the weighted sum of the internal energy (the original cost 
function) and the entropy functional; T and 6 are penalty functionals to en- 
force the constraints of the problem. According to the variational principle of 
Statistical Mechanics, the best approximation to the Gibbs distribution is the 
minimizer of the free energy [18]. Since — TS is a convex functional, the free 
energy is convex at high temperature and the global minimum can be easily 
attained. The solution can then be continuated as temperature is lowered, so 
as to reach a minimum of U. This procedure has shown to be less sensitive to 
local minima than conventional descent methods, and gives results close to the 
ones from Simulated Annealing, while requiring less computational time [19]. 

The equations for the minimum of the free energy are usually called "mean- 
field equations": 

dF dF 

0; T—^- r = (18) 



dp x (x,y,a) dp y (x,y,a) 



After simple calculations, the solution of eqs. (18) is found to have the following 
form: 



exp{-/5 
p x {x, y, a) = — 

E Q '=-L exp 



dU 



dp x (x,y,a) dp x ( 



+ 



8U 



p y (x,y,a) 



exp | f3 | </,),,(.(■. n 



dp x (x,y,a') dp x {x,y,a')\ J 

dU , dT 1 \ 

a)J J 



+ 



Ea'=-L ex P { P [dp y {x U y,a>) + dpy^y,**')] } 



(19) 
(20) 



where the {V} multipliers have been fixed to normalize the distributions, and 
P = 1/T is the inverse temperature. Now we observe that, for each site (x,y) 
on the grid: 



dU 



dU dm x dU dQ x 

+ - - — = Ct- 



dU 



+ a z 



dU 



dp x (a) dm x dp x (a) dQ x dp x (a) dm x (x,y) dQ x (x,y)' 



(21) 
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Analogously, one can easily find: 



dU dU 2 dU 

= ^ 7Z 3 + ^ 7Z ^ ( 22 ) 



dp y (a) dm y (x,y) dQ y (x,yY 
dT 

-^j^ = a[\(x,y) - \(x,y -I)], (23) 

dr 

—— = a [-\{x, y) + X(x - 1, y)\ . (24) 
op y \OL) 

The derivatives of U with respect to the {m} and {Q} variables are reported 
in Appendix A. From these expressions it is clear that the present formulation 
of PU is equivalent to finding the ground state of a finite-size, spin-L Ising 
model with local constraints, and under a non-uniform magnetic field. 

The consistency constraints are written as equations for the {A} field: 



X(x,y) = \(x,y) -b\m x (x,y) + m y (x + l,y) 

-m x (x,y + l)-my(x,y) + I(x,y)\, (25) 

where b is a small constant. 

Equations (19-20) and (25) are the mean-field equations for PU for arbitrary 
L. As already explained, the MFA technique consists in solving iteratively the 
mean-field equations at high temperature (low (3), and then track the solution 
as the temperature is lowered (f3 grows). 

The algorithm can be summarized as follows. The initial distributions give 
the same probability to each value of the correction field, i.e. p x (x,y,a) = 
p y (x,y,a) = tli e inverse temperature is set to Pmin (Pmax is the lowest 
temperature). Then: 

(1) Evaluate {m} and {Q} fields by Eqs. (13-14); 

(2) Iterate Eqs. (19-20); 

(3) Iterate Eq. (25); 

(4) If Eqs. (19-20) or Eq. (25) are not satisfied, goto step 1; 

(5) If (3 < /3max, increase f3 and goto step I. 

The output of this algorithm is a field {itiout} which approximates the average 
of {k} over the global minima of the cost functional R. 

We remark that the output of the algorithm described in [15] satisfies m G 
[—1,1] for each component of {hiout}, whereas the present algorithm satisfies 
the weaker constraint m G [— L, L] and therefore can be used also in the case 
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of high degree of undersampling. The estimate for the true phase gradient is 
(V/)est = A + 27rm uT; the phase pattern / can then be reconstructed by 
(V/) es t as described in [15]. 



4 Experiments 



In this section we describe some experiments we performed to test the effec- 
tiveness of the proposed algorithm. 



In fig. l-(a) a synthetic phase pattern is shown, while in fig. l-(b) the wrapped 
phase pattern is depicted. This test phase pattern has been constructed by 
the following formula: 



= 120 exp [-\r 2 (x, y)(m + /j 2 c(x, y))] , 1 < x, y < 128 



(26) 



with: 



y/(x - 35.5) 2 + (y - 65.5) 2 , 
x — 35.5 

r 

0.01, 
0.0004, 

where <E> is in radians. 

By construction, $ is undersampled: in fig. 2- (a) black pixels represent loca- 
tions where the true correction field k is such that max{k x , k y } = 2, while gray 
pixels represent locations where max{k x , k y } = 1. The residue map is depicted 
in fig. 2-(b). 

We applied the proposed algorithm to unwrap this test pattern. We used 
L = 2, b = 0.05 and the annealing schedule was established as consisting of 
25 temperature values, equally spaced in the interval [/?min = 0.05, Pmax = 
1.5]. The output phase pattern is depicted in fig. 3. The computational time 
was comparable to that corresponding to [15]. The input phase surface was 
perfectly reconstructed. Let us now compare the performance of the MFA 
algorithm with that from LMS [17]. In fig. 4 we show the output of LMS 
applied to the surface of fig. l-(a). Due to severe undersampling, the LMS 
performance is poor. 

We also investigated the robustness of the MFA algorithm with respect to 
noise. In fig. 5- (a) the phase pattern obtained by adding unit- variance Gaus- 
sian noise to the surface of fig. l-(a) is shown. The wrapped phase pattern 



r(x,y) = 
c(x,y) = 

Af2 = 
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is depicted in fig. 5-(b), while in fig. 5-(c) the inconsistencies are shown. The 
output of the MFA algorithm is shown in fig. 6- (a), while in fig. 6-(b) we show 
the phase pattern obtained by re- wrapping the MFA output. The smoothing 
capability of the proposed algorithm appears clearly by comparing figs. 5-(b) 
and 6-(b). 



5 Conclusions 



In this paper we have generalized a previously presented MFA algorithm to 
unwrap phase patterns. This problem is formulated as that of finding the 
ground state of a locally-constrained, spin-L Ising model under a non-uniform 
external field. The present generalization allows processing of noisy and highly 
undersampled input phase fields. The effectiveness of this statistical approach 
to PU has been demonstrated on simulated phase surfaces. Further work will 
be devoted to the estimation of the optimal value of L from the observed 
wrapped phase data. 
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(a) (b) 

Fig. 1. (a) Synthetic phase surface generated via eq. (26); (b) wrapped phase pattern: 
principal phase values span the interval from — ir (black pixels) to +ir (white pixels). 



(a) (b) 

Fig. 2. (a) Degree of undersampling of the phase field depicted in fig. 1: gray pix- 
els represent locations where the absolute phase gradient differs from its estimate 
(wrapping of the principal phase gradient) by one 27r-cycle, black pixels represent 
locations where the difference is 2 cycles; (b) residue map: white pixels are positive 
residues (1 = 1), black pixels are negative residues (I = —1), gray pixels correspond 
to irrotational locations (1 = 0). 
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Fig. 3. Unwrapped phase reconstructed by the proposed algorithm from the wrapped 
phase field shown in fig. l-(b). 



50 r 



40 '- 




Fig. 4. Unwrapped phase reconstructed by the LMS algorithm from the wrapped 
phase field shown in fig. l-(b). 
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Fig. 5. (a) Same synthetic surface as in fig. 1, with added Gaussian noise with unit 
variance; (b) principal phase pattern; (c) residue map. 



1 40 r 




(a) (b) 

Fig. 6. (a) Unwrapped surface reconstructed by the proposed algorithm from the 
wrapped phase field shown in fig. 5-(b); (b) principal phase pattern obtained by 
re- wrapping the solution depicted in (a). 
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A Appendix: Derivatives of the internal energy 



We report here the expressions of the derivatives of the internal energy in the 
Mean-field approximation. 

One easily finds that: 



dU 



dQ x (x,y) 
dU 



= 4 



OQ y (x,y) 



= 4, 



(A.1) 



dU 
dm x (x,y) 



-2m x (x - 1, y) + - [A x (x, y) - A x (x - 1, y)\ + 



7i 



-2m x (x + 1,1/) + - [A x (x, y) - A x (x + 1, y)] + 

71 

-2m x (x, y-l) + - [A x (x, y) - A x (x, y-l)] + 

71 

-2m x (x, y + l) + - [A x (x, y) - A x (x, y + 1)]. 

7T 



(A.2) 



dU 



dm y (x,y) 



= -2m y (x, y - 1) + - [A y (x, y) - A y (x, y-l)] + 

7T 

-2m y (x, y + 1) + - [A y (x, y) - A y (x, y + 1)] + 

7T 

-2m y (x + l,y) + - [A y (x, y) - A y (x + 1, y)\ + 

7T 

-Hi 1 - !> y) + - [A/^' s/) ~ A( x - !> y)l • 

7T 



(A.3) 
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1 (a) Synthetic phase surface generated via eq. (26); (b) wrapped 


phase pattern: principal phase values span the interval from 




—7i (black pixels) to +7r (white pixels). 







2 (a) Degree of undersampling of the phase held depicted in 



fig. 1: gray pixels represent locations where the absolute 




phase gradient differs from its estimate (wrapping of the 




principal phase gradient) by one 27r-cycle, black pixels 


represent locations where the difference is 2 cycles; (b) residue 


map: white pixels are positive residues (/ = 1), black pixels 


are negative residues [1 = —1), gray pixels correspond tc 


> 


irrotational locations (/ = 0). 







3 Unwrapped phase reconstructed by the proposed algorithm 


from the wrapped phase held shown in fig. l-(b). 





1 Unwrapped phase reconstructed by the LMS algorithm from 


the wrapped phase held shown in hg. l-(b). 







5 (a) Same synthetic surface as in hg. 1, with added Gaussian 


noise with unit variance; (b) principal phase pattern; (c) 






residue map. 







6 (a) Unwrapped surlace reconstructed by the proposed 




algorithm from the wrapped phase held shown in hg. 5-(b); (b) 


principal phase pattern obtained by re-wrapping the solution 


depicted in (a). 
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